Final Project - Behind The Scenes

Introduction

Computer science is a field that is growing rapidly in the United States and around the world today. Industry is constantly releasing advancements in computer science and technology is becoming more ingrained into our daily lives with each passing day. Therefore, with an increase in technology usage, the demand for computer scientists has increased in popularity. Due to this fast-growing field, educational institutions and systems are increasing the amount of computer science courses offered in order to train more future computer scientists.

These increases started at the college level, where majoring in computer science is now becoming a widely available option. At Macalester College, it is one of the largest departments for both students and faculty. While the availability of courses at the college level is a good start, there is a big push to have computer science courses offered in K-12 education. Offering computer science courses in elementary and secondary schools provides an opportunity for kids to expose themselves to coding, potentially leading to younger students discovering new interests and actively engaging with computer science earlier. Oftentimes, being exposed to computer science at a younger age can make students more comfortable with the material and the field later on, which can contribute to a more empowered and diverse set of students entering the workforce or higher education. Given the importance of having computer science courses available in K-12 education, we decided to investigate the availability of computer science courses in K-12 school districts in Minnesota. In this project, we explore the connection between a variety of data sets related to this topic, including K-12 computer science course availability in Minnesota, demographic information from the U.S. census, ACT scores, and school district financial information.

To see our final product, visit: https://advanced-data-science-final-project.netlify.app

About the Data Sets

Minnesota Common Course Catalogue

This data set is provided by the Minnesota State Department of Education and is available to the public. The Minnesota Common Course Catalogue (MCCC) presents information about all the courses that are offered across school districts in Minnesota. For our purposes, we are looking at subject area 10, computer and information sciences (K-12), and all categories and classifications. For each course classification, there is information about the district where it is taught, the region of the state, the local course title, whether it is an AP or IB course, if it fulfills a graduation requirement, among others. We joined this data set with information about what category the course classification falls into, such as A - computer literacy, B - management information systems, or C - network systems. The final data set has 1,602 computer and information classes taught across the state in K-12 education.

Link: https://public.education.mn.gov/MDEAnalytics/DataTopic.jsp?TOPICID=84

American Community Survey

We accessed the 2019 American Community Survey (ACS) from the tidycensus package in R. The ACS is an annual demographic survey from the U.S. Census Bureau. We used this data set to find by-district information about race, socioeconomic status, living environments, and other relevant variables for our analyses and visualizations. We hypothesized that computer science course availability would be connected to one or more of these variables, due to the trend where more computer science courses are offered in wealthier areas across the country that have more resources. We are using this data set to further investigate that possibility.

ACT Data

We retrieved this data from the Minnesota State Department of Education, which gives us average composite ACT scores for each school district in Minnesota for a range of years. It has been proven that higher ACT scores are connected to wealth, which would lead us to expect that districts with more resources would have higher average ACT scores. As aforementioned, we predict that computer science course availability is greater in areas with more wealth and resources, therefore we thought these scores may be another interesting variable in predicting or visualizing availability of computer science courses by district. In the same way we expect certain demographic information to be correlated with computer science availability, we are hoping to use average ACT score information to present another method of showing the connection.

Link: https://public.education.mn.gov/MDEAnalytics/DataTopic.jsp?TOPICID=87

Annual Survey of School System Finances

Another important set of data that we used from the U.S. Census Bureau is the Annual Survey of School System Finances (ASSSF). This data contains information about 2018 financial activity of public elementary and secondary school systems for all states across the country. In particular, we were interested in discovering if funding, revenue, or spending for a district was correlated to availability of computer science courses. Due to the fact that public school funding often comes from income taxes, wealthier areas provide more funding for their local schools.

Link: https://www.census.gov/programs-surveys/school-finances.html

Minnesota School Districts

This data set includes the name of all school districts in Minnesota as well as their district numbers, which was created in 2015. This data came from the Minnesota State Department of Human Services. We joined the district numbers from this data set to the American Community Survey data, which only had the district names. Accessing district numbers and names was integral in connecting our various data sets because the common piece of information between our ideas is district information.

Link: https://www.dhs.state.mn.us/main/groups/county_access/documents/pub/dhs16_193591.pdf

Important Packages and Data Sets

Load Packages

library(tidycensus) # for getting the census data
library(tidyverse) # for data cleaning and visualization
library(sf) # for mapping
library(tidymodels)
library(naniar)
library(DALEX)
library(DALEXtra)
library(vip)
library(gridExtra)
library(plotly)

Load Data Sets

act <- read_csv("act_data.csv") # ACT
districts <- read_csv("districts.csv") # MN School Districts
course_catalog <- read_csv("course_catalog.csv") # MCCC
ss <- readxl::read_excel("elsec18t.xls") # ASSSF

Data Cleaning

MCCC Data

The data from the Department of Education was neatly split into spreadsheets based on course classifications. In order to use this data more easily, we manually combined all of the spreadsheets into one data set with all of the courses offered. We had to manually add in the columns for course classification and category in order to create one data set. Furthermore, in order to connect to the other data sets, we needed to have the district number as its own variable. In this data set, the district name and number were combined as one variable, so we used regular expressions to split the name and number. This was fairly simple, due to the district variable consistently having the number as the last 7 digits.

course_catalog_predictors <-
course_catalog %>%
  mutate(DistNum = as.integer(str_sub(`District`, -7,-4))) %>%
  group_by(DistNum) %>%
  summarise(TotalClasses = sum(n()),
            NumCat = length(unique(Category)))

ACS Data

The ACS data was already in a tidy format, which was convenient for our purposes. We focused on narrowing down variables to the ones that were important to our project and renaming them so that we could easily understand the meaning of the variable. Many of the districts had the word ‘Minnesota’ in it, so we removed that because all of our data is about Minnesota. We also created some percentage variables because they would be more meaningful for our exploration than the raw numbers that tidycensus provided us with.

# census_api_key("key")

mn_2019_census <- get_acs(state = "MN", 
                         geography = "school district (unified)",
                         variables = c("B01003_001",
                                       "B01001_001",
                                       "B01001_004",
                                       "B01001_005",
                                       "B01001_006",
                                       "B02001_001",
                                       "B02001_002",
                                       "B02001_003",
                                       "B02001_005",
                                       "B02001_006",
                                       "B02001_008",
                                       "B19013_001",
                                       "B09010_001",
                                       "B09010_002", 
                                       "B25031_001",
                                       "B25027_001",
                                       "B25027_002",
                                       "B28002_001",
                                       "B28002_002"),
                         geometry = TRUE,
                         year = 2019)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================| 100%

mn_2019_census_nogeom <- mn_2019_census %>%
  select(-`moe`) %>% #get rid of moe variable(not sure what it is)
  spread(variable, estimate) %>% #give each variable a column
  rename(total_pop = B01003_001,
         male_5to9 = B01001_004,
         male_10to14 = B01001_005,
         male_15to17 = B01001_006,
         race_white_only = B02001_002,
         race_Black_only = B02001_003,
         race_Asian_only = B02001_005,
         race_pacific_islander_only = B02001_006,
         at_least_two_races = B02001_008,
         med_household_inc_12mo = B19013_001,
         SSI_pubassist_foodstamps = B09010_002,
         med_gross_rent = B25031_001,
         house_units_w_mortgage = B25027_002,
         internet_subscrip_in_house = B28002_002
         ) %>% #rename variables
  mutate(District = str_extract(`NAME`, "[^,]+")) %>% #delete ", Minnesota" after district name
  select(-NAME) %>% #remove `NAME` column
  mutate(perc_male_5to9 = male_5to9/B01001_001,
         perc_male_10to14 = male_10to14/B01001_001,
         perc_male_15to17 = male_15to17/B01001_001,
         perc_white_only = race_white_only/B02001_001,
         perc_black_only = race_Black_only/B02001_001,
         perc_asian_only = race_Asian_only/B02001_001,
         perc_pacific_islander_only = race_pacific_islander_only/B02001_001,
         perc_at_least_two_races = at_least_two_races/B02001_001,
         perc_SSI_pubassist_foodstamps = SSI_pubassist_foodstamps/B09010_001,
         perc_house_units_w_mortgage = house_units_w_mortgage/B25027_001,
         perc_internet_subscription = internet_subscrip_in_house/B28002_001
         ) %>% #create percentage variables
  select(-B01001_001,-B02001_001, -B09010_001, -B25027_001, -B28002_001 ) #deselect variables used to calculate % but aren't important for our visualizations

ACT Data

This data set included a range of years, so we filtered the graduation year to our year of interest: 2018. This data set also originally had extra variables that were irrelevant to our interests, so our cleaning process involved removing predictors so that we were only left with district name and average composite ACT score. The district name included the district number, so in order to create two separate values we used regular expressions. The district number was consistently located at the same spot in the name, so splitting them was relatively straightforward.

act_clean <-
act %>%
  filter(`Grad Year` == 2018,
         `Analysis Level` == "District",
         `District Name` != "MINNESOTA DEPT OF EDUCATION",
         `Avg Comp` != ".") %>%
  select(`District Name`, `Avg Comp`) %>%
  mutate(`Dist Num` = as.integer(str_extract_all(`District Name`, "[:digit:]+"))) %>%
  mutate(`Dist Num`= replace(`Dist Num`, `District Name`=="GRANADA-HUNTLEY-EAST CHAIN SD", 2536),
         `Dist Num`= replace(`Dist Num`, `District Name`=="ALBERT LEA AREA SCHOOLS", 241),
         `Dist Num`= replace(`Dist Num`, `District Name`=="ATWATER-COSMOS-GROVE CITY SD", 2396),
         `Dist Num`= replace(`Dist Num`, `District Name`=="BOLD SCHOOL DISTRICT", 2534),
         `Dist Num`= replace(`Dist Num`, `District Name`=="BUFFALO HANOVER MONTROSE SD", 0877),
         `Dist Num`= replace(`Dist Num`, `District Name`=="BUFFALO LAKE-HECTOR-STEWART SD", 2159),
         `Dist Num`= replace(`Dist Num`, `District Name`=="EASTERN CARVER CO SCHOOLS", 112),
         `Dist Num`= replace(`Dist Num`, `District Name`=="MANKATO AREA PUBLIC SCHS", 77),
         `Dist Num`= replace(`Dist Num`, `District Name`=="MINNEAPOLIS PUBLIC SCH DIST", 9991),
         `Dist Num`= replace(`Dist Num`, `District Name`=="MINNEOTA PUBLIC SCHOOLS", 414),
         `Dist Num`= replace(`Dist Num`, `District Name`=="PLAINVIEW ELGIN MILLVILLE SD", 2899),
         `Dist Num`= replace(`Dist Num`, `District Name`=="ROYALTON PUBLIC SCHOOLS", 485),
         `Dist Num`= replace(`Dist Num`, `District Name`=="WASECA PUBLIC SCHOOL DIST", 829),
         `Dist Num`= replace(`Dist Num`, `District Name`=="WAUBUN-OGEMA-WHITE EARTH PSD", 435),
         `Dist Num`= replace(`Dist Num`, `District Name`=="WINONA AREA PUBLIC SCHOOL DIST", 861))

Annual Survey of School System Finances

Fortunately, we did not have to clean this data set. Originally, this data set included information about all of the states, but the process of data joining narrowed down the data to just contain information about Minnesota. Below is the code that allowed us to view the data entries related to Minnesota before going through the process of data joining.

ss %>%
  # filter for MN
  filter(str_detect(NCESID, "^27")) %>%
  arrange(NCESID) %>%
  head(10) %>%
  select(NCESID)

Minnesota School Districts

This data set originally had four categories: district name, district number, start date, and end date. To prepare this data set for joining, we removed the start and end date variables. In order to join it with the ACS data set, we had to compare the district names in both to ensure that they aligned. Only 16 district names didn’t match, which was a small number, so we manually changed the district names to match the names in the ACS data set.

Data Joining

Part of the work we did after finding, importing, and cleaning our various data sets was joining them together. An important piece for this was ensuring that the district names and numbers aligned, and most of that was completed in the previous step. Here is our code for joining the data sets:

# read in data sets
## annual survey of school system finances (2018, national)
ss <- readxl::read_excel("elsec18t.xls")

## selected variables from the ACS (2019, MN)
mn_acs <- st_read("mn_2019_census/mn_2019_census.shp", 
                  geometry_column = "geometry", 
                  fid_column_name = "geometry")

## Reading layer `mn_2019_census' from data source `/Users/anaelkuperwajs/STAT 494/STAT494-Final-Project/mn_2019_census/mn_2019_census.shp' using driver `ESRI Shapefile'
## Simple feature collection with 324 features and 28 fields (with 1 geometry empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -97.23921 ymin: 43.49936 xmax: -89.49174 ymax: 49.38436
## Geodetic CRS:  NAD83

mn_acs <- as.data.frame(mn_acs)
# rename the columns
names(mn_acs) <- c("GEOID", "male_5to9", "male_10to14", "male_15to17",
                   "total_pop", "race_white_only", "race_Black_only", "race_Asian_only",
                   "race_pacific_islander_only", "at_least_two_races", "SSI_pubassist_foodstamps",
                   "med_household_inc_12mo", "house_units_w_mortgage", "med_gross_rent",
                   "internet_subscrip_in_house", "District", "perc_male_5to9", "perc_male_10to14",
                   "perc_male_15to17", "perc_white_only", "perc_black_only", "perc_asian_only",
                   "perc_pacific_islander_only", "perc_at_least_two_races", "perc_SSI_pubassist_foodstamps",
                   "perc_house_units_w_mortgage", "perc_internet_subscription",
                   "District.Nbr", "geometry")
# move the district name and number cols to the front
mn_acs <- mn_acs %>%
  relocate(District, District.Nbr, .after = GEOID)

## ACT scores (MN)
mn_act <- read.csv("act_clean.csv")[-1] # remove the redundant index column created when reading by csv

Join by NCES ID for school districts.

# look at this column from the school survey
ss %>%
  # filter for MN 
  filter(str_detect(NCESID, "^27")) %>%
  arrange(NCESID) %>%
  head(10) %>%
  select(NCESID)

# and from the ACS
mn_acs %>%
  arrange(GEOID) %>%
  head(10) %>%
  select(GEOID)

# left join on ACS data
mn_acs_ss <- mn_acs %>%
  left_join(ss, by = c("GEOID" = "NCESID"))

Let’s look at the name columns from the two data sets to verify that they point to the same district.

mn_acs_ss %>%
  select(District, NAME)

Now let’s see which rows from the ACS and school survey didn’t get a match.

# ACS
mn_acs_ss %>%
  filter(is.na(NAME))

# school survey
ss %>%
  filter(str_detect(NCESID, "^27")) %>% # filter for MN
  anti_join(mn_acs, by = c("NCESID" = "GEOID"))

So only the row “Remainder of Minnesota” from the ACS didn’t have a match, which is reasonable. There are 67 school districts from the school survey that did not join to the ACS, but perhaps they fall into the “Remainder of Minnesota” category. A lot of these districts also had 0 for enrollment.

To join this data set with the ACT data, we’ll have to use the school district numbers as the ACT data does not contain the NCES ID.

# look at this variable from each data set to see their format
mn_acs_ss %>%
  arrange(District.Nbr) %>%
  select(District.Nbr)

mn_act %>%
  arrange(Dist.Num) %>%
  select(Dist.Num)

It appears that we’re missing ACT scores for 11 school districts. Nevertheless, the formats do match so let’s join the data sets.

mn_acs_ss_act <- mn_acs_ss %>%
  mutate(District.Nbr = as.integer(District.Nbr)) %>%
  left_join(mn_act, by = c("District.Nbr" = "Dist.Num"))

Finally, before we can export the data, we will remove the redundant names and ID columns from the school survey and ACT data sets as we’ll stick to their format from the ACS.

mn_acs_ss_act <- mn_acs_ss_act %>%
  select(-c("IDCENSUS", "NAME", "District.Name"))

Export the data:

#write.csv(mn_acs_ss_act, "mn_acs_ss_act.csv")

st_write(mn_acs_ss_act, "mn_acs_ss_act/mn_acs_ss_act.shp", 
         append=FALSE, 
         fid_column_name = "geometry")

## Deleting layer `mn_acs_ss_act' using driver `ESRI Shapefile'
## Writing layer `mn_acs_ss_act' to data source `mn_acs_ss_act/mn_acs_ss_act.shp' using driver `ESRI Shapefile'
## Writing 324 features with 92 fields and geometry type Multi Polygon.

Join with predictor data, replace all missing data from the course catalog with 0, and remove the redundant variable. This is the final step!

mn_acs_ss_act_pred <- mn_acs_ss_act %>%
  # mutate(District.Nbr = as.integer(District.Nbr)) %>% ## no longer necessary since they're already ints
  ## changed to left join to ignore rows without matches from the course catalog
  left_join(course_catalog_predictors, by = c("District.Nbr" = "DistNum")) %>% 
  mutate(TotalClasses = ifelse(is.na(TotalClasses), 0, TotalClasses),
         NumCat = ifelse(is.na(NumCat), 0, NumCat)) 

Mapping and Shape Files

While we had all the data necessary to create useful graphs and tables, one of our goals was to create visualizations and maps. Our data is focused on the state of Minnesota and all of it is segmented by school districts. Therefore, we believed that an interesting and informative way to present this data is through maps. The piece that was missing for us to accomplish this was the shape files for the districts. Shape files are notoriously tricky, and implementing this was a challenge.

Once we added a column for the geometry of each district (which was necessary for mapping), we could no longer use the write and read_csv functions that we were so accustomed to. Unfortunately, we found this out the hard way. When we exported our final data set using write_csv() and read it in another file with read_csv(), our data values were completely altered. We hypothesized that this could be due to how each geometry feature is saved as a list of coordinates separated by commas, which would interfere with the csv (comma-separated values) format. Fortunately, we discovered that st_write() and st_read() from the sf package were the functions we needed to use instead. These were specifically made for exporting and importing files with simple geometry features. One issue we encountered with exporting shape files is that certain column names were abbreviated in the process. For example, male_5to9 became mal_5t9 and total_pop became totl_pp. In the end, we decided that the simplest and fastest solution was to rename our variables back to their original names after we read in the data.

##### RUN THIS CHUNK WHENEVER YOU READ IN THE DATA
mn_acs_ss_act_pred <- st_read("mn_acs_ss_act_pred/mn_acs_ss_act_pred.shp",
                geometry_column = "geometry",
                fid_column_name = "geometry",
                quiet = TRUE)
# rename the census variables because their names were reformatted...
names(mn_acs_ss_act_pred)[1:28] <- 
  c("GEOID", "District", "District.Nbr", "male_5to9", "male_10to14", "male_15to17",
    "total_pop", "race_white_only", "race_Black_only", "race_Asian_only",
    "race_pacific_islander_only", "at_least_two_races", "SSI_pubassist_foodstamps",
    "med_household_inc_12mo", "house_units_w_mortgage", "med_gross_rent",
    "internet_subscrip_in_house", "perc_male_5to9", "perc_male_10to14",
    "perc_male_15to17", "perc_white_only", "perc_black_only", "perc_asian_only",
    "perc_pacific_islander_only", "perc_at_least_two_races", "perc_SSI_pubassist_foodstamps",
    "perc_house_units_w_mortgage", "perc_internet_subscription")
names(mn_acs_ss_act_pred)[93] <- c("TotalClasses")

Visualizations and Plotly

Although we knew our maps up to this point were telling, we wanted to add another element to our project. We decided to make the maps interactive and chose to use plotly to achieve this. The plotly functions had an issue with our final data set and would not map unless we removed all rows with NA; converting the NA values to 0 did not fix the issue, and neither did converting the ‘N’ values we found in two variable columns to 0. In total, dropping the NA rows removed 11 districts. With more time, this issue should be explored further. Something that is important to note is that the districts that were not mapped still show interactivity based on whichever district is closest to the pointer.

Below we have two examples from our final product: one of a static graph and one of a plotly graph. Most of our visualizations and maps were built in the same way, so we did not include all of the visualization code in this document.

Static Graph

#population
ggplot(mn_acs_ss_act_pred, aes(fill = total_pop)) + 
  geom_sf(color = NA) + 
  scale_fill_viridis_c(option = "magma", labels = scales::comma) +
  labs(title = "Total Population", fill = "") +
  theme(axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.background = element_rect(fill = "white"),
        plot.title = element_text(hjust = 0.5))

Plotly

mn_acs_ss_act_pred_viz <- 
  mn_acs_ss_act_pred %>%  drop_na()

ggplotly(
  ggplot(mn_acs_ss_act_pred_viz) +
  geom_sf(aes(fill = TotalClasses, text = paste(District)),  size = .1, color = "white") +
  scale_fill_viridis_c(option = "B") +
  labs(title = "Total CS Courses Offered", fill = "") +
  theme(axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.background = element_rect(fill = "white"),
        plot.title = element_text(hjust = 0.5),
        title =element_text(size=12, face='bold'))
) 

Modeling

The goal of our model was to predict the number of computer science courses offered per district in Minnesota based on the variables from our ACS, ASSSF, and ACT data.

We first read in the data and cleaned up variable names, as they were abbreviated after being exported to shape file format.

# read in the data
mn <- st_read("mn_acs_ss_act_pred/mn_acs_ss_act_pred.shp",
              geometry_column = "geometry",
              fid_column_name = "geometry")

## Reading layer `mn_acs_ss_act_pred' from data source `/Users/anaelkuperwajs/STAT 494/STAT494-Final-Project/mn_acs_ss_act_pred/mn_acs_ss_act_pred.shp' using driver `ESRI Shapefile'
## Simple feature collection with 324 features and 94 fields (with 1 geometry empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -97.23921 ymin: 43.49936 xmax: -89.49174 ymax: 49.38436
## Geodetic CRS:  NAD83

# rename the census variables because their names were reformatted...
names(mn)[1:28] <- 
  c("GEOID", "District", "District.Nbr", "male_5to9", "male_10to14", "male_15to17",
    "total_pop", "race_white_only", "race_Black_only", "race_Asian_only",
    "race_pacific_islander_only", "at_least_two_races", "SSI_pubassist_foodstamps",
    "med_household_inc_12mo", "house_units_w_mortgage", "med_gross_rent",
    "internet_subscrip_in_house", "perc_male_5to9", "perc_male_10to14",
    "perc_male_15to17", "perc_white_only", "perc_black_only", "perc_asian_only",
    "perc_pacific_islander_only", "perc_at_least_two_races", "perc_SSI_pubassist_foodstamps",
    "perc_house_units_w_mortgage", "perc_internet_subscription")
names(mn)[93] <- c("TotalClasses")

Data Exploration

# look at the variables and their types
str(mn)

## Classes 'sf' and 'data.frame':   324 obs. of  95 variables:
##  $ GEOID                        : chr  "2700001" "2700005" "2700006" "2700007" ...
##  $ District                     : chr  "Mountain Iron-Buhl School District" "United South Central School District" "Maple River School District" "Kingsland Public School District" ...
##  $ District.Nbr                 : int  712 2134 2135 2137 2142 2143 2144 2149 2154 2155 ...
##  $ male_5to9                    : num  146 220 207 133 370 176 625 344 241 266 ...
##  $ male_10to14                  : num  104 160 180 176 493 201 991 272 252 237 ...
##  $ male_15to17                  : num  91 126 133 91 239 93 530 179 132 144 ...
##  $ total_pop                    : num  4545 5790 5991 5023 17435 ...
##  $ race_white_only              : num  4251 5517 5802 4773 15739 ...
##  $ race_Black_only              : num  63 10 26 7 71 21 40 32 30 68 ...
##  $ race_Asian_only              : num  81 43 17 66 49 12 230 54 12 0 ...
##  $ race_pacific_islander_only   : num  0 0 0 0 19 3 0 0 0 0 ...
##  $ at_least_two_races           : num  104 154 137 133 604 70 369 126 111 183 ...
##  $ SSI_pubassist_foodstamps     : num  238 313 196 139 525 151 421 224 239 510 ...
##  $ med_household_inc_12mo       : num  50550 53699 64451 59087 57221 ...
##  $ house_units_w_mortgage       : num  756 1009 1074 1031 3604 ...
##  $ med_gross_rent               : num  625 640 636 675 709 769 929 616 677 670 ...
##  $ internet_subscrip_in_house   : num  1593 1891 1994 1689 6243 ...
##  $ perc_male_5to9               : num  0.0321 0.038 0.0346 0.0265 0.0212 ...
##  $ perc_male_10to14             : num  0.0229 0.0276 0.03 0.035 0.0283 ...
##  $ perc_male_15to17             : num  0.02 0.0218 0.0222 0.0181 0.0137 ...
##  $ perc_white_only              : num  0.935 0.953 0.968 0.95 0.903 ...
##  $ perc_black_only              : num  0.01386 0.00173 0.00434 0.00139 0.00407 ...
##  $ perc_asian_only              : num  0.01782 0.00743 0.00284 0.01314 0.00281 ...
##  $ perc_pacific_islander_only   : num  0 0 0 0 0.00109 ...
##  $ perc_at_least_two_races      : num  0.0229 0.0266 0.0229 0.0265 0.0346 ...
##  $ perc_SSI_pubassist_foodstamps: num  0.285 0.222 0.137 0.124 0.184 ...
##  $ perc_house_units_w_mortgage  : num  0.505 0.49 0.552 0.572 0.513 ...
##  $ perc_internet_subscription   : num  0.756 0.771 0.823 0.78 0.779 ...
##  $ CONUM                        : chr  "27137" "27043" "27013" "27045" ...
##  $ CSA                          : chr  "N" "N" "359" "462" ...
##  $ CBSA                         : chr  "20260" "N" "31860" "40340" ...
##  $ ENROLL                       : num  507 707 927 557 2007 ...
##  $ TOTALRE                      : num  8146 12242 13103 8078 39951 ...
##  $ TFEDREV                      : num  442 554 489 374 1860 503 835 848 545 955 ...
##  $ FEDRCOM                      : num  175 170 146 160 506 83 197 194 168 321 ...
##  $ FEDRSPE                      : num  0 0 0 0 0 228 0 0 0 0 ...
##  $ FEDRNUT                      : num  133 215 213 135 575 166 474 283 186 483 ...
##  $ FEDROTH                      : num  134 169 130 79 779 26 164 371 191 151 ...
##  $ TSTREV                       : num  5171 7246 9420 5269 22890 ...
##  $ STRFORM                      : num  3618 5554 7178 4295 17628 ...
##  $ STRSPEC                      : num  253 818 1379 339 2493 ...
##  $ STRTRAN                      : num  0 7 4 0 0 17 929 27 0 3 ...
##  $ STROTHR                      : num  1300 867 859 635 2769 ...
##  $ TLOCREV                      : num  2533 4442 3194 2435 15201 ...
##  $ LOCRTAX                      : num  294 3509 1903 1685 6434 ...
##  $ LOCRPRO                      : num  294 3509 1903 1685 6434 ...
##  $ LOCRPAR                      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ LOCRCIC                      : num  1642 45 63 2 6756 ...
##  $ LOCROSC                      : num  37 4 151 48 406 103 1 32 54 117 ...
##  $ LOCRCHA                      : num  132 383 855 510 643 ...
##  $ LOCROTH                      : num  428 501 222 190 962 ...
##  $ TOTALEX                      : num  26135 10618 12460 7379 37999 ...
##  $ TCURSPN                      : num  6871 9142 11697 6468 29197 ...
##  $ TSALWAG                      : num  4021 5639 6706 3799 16988 ...
##  $ TEMPBEN                      : num  1496 1927 1789 977 6655 ...
##  $ TCURINS                      : num  3915 5729 6241 3755 18391 ...
##  $ TCURISA                      : num  2601 3974 4431 2515 11543 ...
##  $ TCURIBE                      : num  954 1221 1159 634 4455 ...
##  $ TCURSSV                      : num  2547 2722 4290 2054 9060 ...
##  $ TCURSPU                      : num  128 206 428 107 495 ...
##  $ TCURSST                      : num  167 237 533 229 813 ...
##  $ TCURSGE                      : num  369 523 239 208 1075 ...
##  $ TCURSSC                      : num  352 446 523 332 1454 ...
##  $ TCURSOT                      : num  1531 1310 2567 1178 5223 ...
##  $ TCURONO                      : num  409 691 1166 659 1746 ...
##  $ TCAPOUT                      : num  18168 561 473 333 5900 ...
##  $ TPAYOTH                      : num  224 186 286 347 947 ...
##  $ TINTRST                      : num  872 729 4 231 1955 ...
##  $ DEBTOUT                      : num  29340 24889 0 18095 56845 ...
##  $ LONGISS                      : num  0 0 0 0 0 ...
##  $ LONGRET                      : num  417 1239 455 767 3715 ...
##  $ PCTTOTA                      : num  100 100 100 100 100 100 100 100 100 100 ...
##  $ PCTFTOT                      : num  5.4 4.5 3.7 4.6 4.7 4.9 1.8 4.1 4.1 6.8 ...
##  $ PCTFCOM                      : num  2.1 1.4 1.1 2 1.3 0.8 0.4 0.9 1.3 2.3 ...
##  $ PCTSTOT                      : num  63.5 59.2 71.9 65.2 57.3 66.7 70.3 65.8 64.3 76.1 ...
##  $ PCTSFOR                      : num  44.4 45.4 54.8 53.2 44.1 52.7 50.2 41.4 51.7 58.6 ...
##  $ PCTLTOT                      : num  31.1 36.3 24.4 30.1 38 28.5 27.9 30.1 31.6 17.1 ...
##  $ PCTLTAX                      : num  3.6 28.7 14.5 20.9 16.1 18.5 18.1 18.5 3.5 9.1 ...
##  $ PCTLOTH                      : num  20.6 0.4 1.6 0.6 17.9 1.5 0.4 0.6 22 1 ...
##  $ PCTLCHA                      : num  1.6 3.1 6.5 6.3 1.6 5.2 6.5 4 2.5 4.8 ...
##  $ PPCSTOT                      : num  13215 12646 12006 11056 14162 ...
##  $ PPSALWG                      : num  7931 7976 7234 6820 8464 ...
##  $ PPEMPBE                      : num  2951 2726 1930 1754 3316 ...
##  $ PPITOTA                      : num  7722 8103 6732 6741 9163 ...
##  $ PPISALW                      : num  5130 5621 4780 4515 5751 ...
##  $ PPIEMBE                      : num  1882 1727 1250 1138 2220 ...
##  $ PPSTOTA                      : num  5024 3850 4628 3688 4514 ...
##  $ PPSPUPI                      : num  252 291 462 192 247 329 318 152 159 215 ...
##  $ PPSSTAF                      : num  329 335 575 411 405 473 724 333 605 440 ...
##  $ PPSGENA                      : num  728 740 258 373 536 290 268 185 457 271 ...
##  $ PPSSCHA                      : num  694 631 564 596 724 698 467 299 424 441 ...
##  $ Avg_Cmp                      : num  19.1 20.5 21.1 21 20.1 19.1 21.8 20.4 19.6 19.7 ...
##  $ TotalClasses                 : num  0 1 6 9 7 6 12 7 2 3 ...
##  $ NumCat                       : num  0 1 3 1 2 1 3 3 2 2 ...
##  $ geometry                     :sfc_MULTIPOLYGON of length 324; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:88, 1:2] -92.8 -92.8 -92.8 -92.8 -92.8 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:94] "GEOID" "Distrct" "Dstrc_N" "mal_5t9" ...

After looking at the variables and their types, we noticed that we were dealing with predominantly numerical variables. GEOID, District, District.Nbr, CONUM, CSA, and CBSA are ID variables, so we will have to specify that in the recipe so that they won’t be included in the model.

We decided to drop the geometry column and convert the data to a regular data frame format for this part because it was not necessary for the models and including it led to issues running some functions.

mn <- mn %>%
  st_drop_geometry()

We plotted the distributions of our variables to see what kind of data transformation we would need to perform.

Below are sets of our predictors, first from the Annual Survey of School System Finances:

mn %>%
  select(-GEOID) %>%
  # variables from the Annual Survey of School System Finances are in all caps
  select(matches("^[A-Z]{4,}", ignore.case = FALSE)) %>%
  select(where(is.numeric)) %>%
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "values") %>%
  ggplot(aes(x = values)) +
  geom_histogram() +
  facet_wrap(vars(variable),
             scales = "free")

Now from the ACS and ACT data:

mn %>%
  # remove the ID variables
  select(-c("GEOID", "District", "District.Nbr")) %>%
  # these variables have at least one lower case character
  select(matches("[a-z]", ignore.case = FALSE)) %>%
  # filter for numerical variables
  select(where(is.numeric)) %>%
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "values") %>%
  ggplot(aes(x = values)) +
  geom_histogram() +
  facet_wrap(vars(variable),
             scales = "free")

A lot of the data are right-skewed, which makes sense as many variables are raw counts.

We also checked for missing data:

mn %>% 
  add_n_miss() %>% 
  count(n_miss_all)

We had one observation missing 65 variables, and this was the “Remainder of Minnesota” observation from the ACS. The NAs caused issues when we tried to fit a LASSO and random forest model, so we removed them beforehand.

mn <- mn %>%
  drop_na()

Model Recipe

From our data exploration, these were the things we had to do in the recipe:

• Log-transform most of the variables (except for the ones that start with PCT from the School Survey as they’re percentages and PP as they’re spending per pupil and not right-skewed).
• Use percentages for ACS variables.
• Normalize all numerical variables.
• Ignore PCTTOTA and LOCRPAR as all variables have the same value.
• Make the ID variables evaluative (i.e. not included in modeling).

Split data

We split the data into a training and testing set and also created cross-validation folds for evaluation.

# split the data into a training and test set
mn_split <- initial_split(mn, prop = .75)
mn_training <- training(mn_split)
mn_testing <- testing(mn_split)

# create cross-validation folds
mn_cv <- vfold_cv(mn_training, v = 5)

Preprocess data

Here is the recipe for our models:

mn_recipe <- 
  recipe(TotalClasses ~ ., data = mn_training) %>%
  # ignore observations with missing data (necessary for LASSO mod)
  step_naomit(everything(), skip = TRUE) %>%
  # remove variables
  step_rm(
    # this one can be considered as response variable itself
    NumCat, 
    # all variables have the same value for these two
    PCTTOTA, LOCRPAR, 
    # raw counts from ACS
    matches("[a-z]", ignore.case = FALSE),
    -starts_with("perc"),
    -total_pop,
    -Avg_Cmp,
    -TotalClasses
    ) %>%
  # log-transform 
  step_log(
    # total population
    total_pop,
    # spending / revenue variables from the school survey
    ## ignore those that start with P since they're percentages / spending per student
    matches("^[A-OQ-Z]{4,}", ignore.case = FALSE), 
    ## ignore ID variables as well
    -GEOID, -CONUM, -CBSA,
    # some variables have 0s which will produce NaNs when log-transformed
    offset = 1) %>% 
  # make ID variables evaluative (not included in modeling)
  update_role(
    all_of(c("GEOID",
             "District",
             "District.Nbr",
             "CONUM",
             "CSA",
             "CBSA")),
    new_role = "evaluative") %>%
  # make integers numeric
  step_mutate_at(is.integer, fn = as.numeric) %>%
  # normalize numerical variables
  step_normalize(all_predictors())

Here is what the data looked like post-transformation:

mn_recipe %>%
  prep(mn_training) %>%
  juice()

Model Fitting

Regular linear regression

We first tested a regular linear regression model and looked at the table of coefficients. Since we had an overwhelming number of predictors, we assumed beforehand that this model would not perform well due to overfitting.

# define the model type
mn_linear_mod <-
  linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")

# set up the workflow
mn_lm_wf <-
  workflow() %>%
  add_recipe(mn_recipe) %>%
  add_model(mn_linear_mod)

# fit the model
mn_lm_fit <-
  mn_lm_wf %>%
  fit(mn_training)

# display the results
mn_lm_fit %>% 
  pull_workflow_fit() %>% 
  tidy() %>% 
  mutate(across(where(is.numeric), ~round(.x,3))) 

LASSO

To deal with the massive number of predictors, we switched to LASSO to shrink coefficients to zero and thereby eliminate insignificant variables from the model.

# define the model type
mn_lasso_mod <-
  linear_reg(mixture = 1) %>%
  set_engine("glmnet") %>%
  set_args(penalty = tune()) %>%
  set_mode("regression")

# set up the workflow
mn_lasso_wf <-
  workflow() %>%
  add_recipe(mn_recipe) %>%
  add_model(mn_lasso_mod)

# set up penalty grid for tuning
penalty_grid <- grid_regular(penalty(),
                             levels = 20)

# tune the parameter
mn_lasso_tune <-
  mn_lasso_wf %>%
  tune_grid(
    resamples = mn_cv,
    grid = penalty_grid
  )

We chose the best parameter based on RMSE and finalized the workflow / model. We then looked at the variables that were retained by LASSO.

# show the best penalty parameter
mn_lasso_tune %>% 
  show_best(metric = "rmse")

# select best parameter by smallest rmse
(best_param <- mn_lasso_tune %>% 
    select_best(metric = "rmse"))

# finalize workflow
mn_lasso_final_wf <- mn_lasso_wf %>% 
  finalize_workflow(best_param)

# fit final model
mn_lasso_final_mod <-
  mn_lasso_final_wf %>%
  fit(data = mn_training)

# look at the table of coefficients
mn_lasso_final_mod %>% 
  pull_workflow_fit() %>% 
  tidy()  %>%
  # filter for predictors with non-zero coefficients 
  filter(estimate != 0)

Random forest

Our last model candidate was a random forest model.

# define the model type
mn_rf_mod <- 
  rand_forest(mtry = 23, # ~1/3 of predictors 
              min_n = 5, 
              trees = 200) %>% 
  set_mode("regression") %>% 
  set_engine("ranger")

# set up the workflow
mn_rf_wf <-
  workflow() %>%
  add_recipe(mn_recipe) %>%
  add_model(mn_rf_mod)

# fit the model
mn_rf_fit <- 
  mn_rf_wf %>% 
  fit(mn_training)

Model Evaluation and Comparison

# fit model with best tuning parameter(s) to training data and apply to test data
mn_lm_test <- 
  mn_lm_wf %>% 
  last_fit(mn_split)
mn_lasso_test <- 
  mn_lasso_final_wf %>% 
  last_fit(mn_split)
mn_rf_test <-
  mn_rf_wf %>%
  last_fit(mn_split)

Of our three models, the random forest performed best, followed by LASSO, and then regular linear regression.

# collect metrics for model applied to test data
mn_lm_test %>%
  collect_metrics

mn_lasso_test %>% 
  collect_metrics()

mn_rf_test %>% 
  collect_metrics()

Residuals

Since the RMSE from the regular lm model was rather high, we decided to only compare the LASSO and random forest models moving forward. We computed their overall performance metrics and looked at the residuals.

lasso_explain <- 
  explain_tidymodels(
    model = mn_lasso_final_mod,
    data = mn_training %>% select(-TotalClasses), 
    y = mn_training %>%  pull(TotalClasses),
    label = "lasso"
  )

## Preparation of a new explainer is initiated
##   -> model label       :  lasso 
##   -> data              :  235  rows  93  cols 
##   -> target variable   :  235  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.2 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  2.983908 , mean =  4.855319 , max =  43.34407  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -18.45726 , mean =  0.000000000000004801915 , max =  18.7106  
##   A new explainer has been created! 

rf_explain <- 
  explain_tidymodels(
    model = mn_rf_fit,
    data = mn_training %>% select(-TotalClasses), 
    y = mn_training %>%  pull(TotalClasses),
    label = "rf"
  )

## Preparation of a new explainer is initiated
##   -> model label       :  rf 
##   -> data              :  235  rows  93  cols 
##   -> target variable   :  235  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.2 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  0.64875 , mean =  4.921705 , max =  43.57892  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -8.43225 , mean =  -0.06638546 , max =  10.42108  
##   A new explainer has been created! 

# get overall performance metrics
lasso_mod_perf <- model_performance(lasso_explain)
rf_mod_perf <-  model_performance(rf_explain)

Here are tables of their performance metrics:

# LASSO
data.frame(lasso_mod_perf$measures)

# random forest
data.frame(rf_mod_perf$measures)

Here is the distribution of the residuals:

plot(lasso_mod_perf,
     rf_mod_perf, 
     geom = "boxplot")

Variable importance

In the end, the random forest showed to greatly outperform the LASSO. Our final step was to look at the variable importance plot from this final model.

set.seed(1) 
# create explainer
rf_explain <- 
  explain_tidymodels(
    model = mn_rf_fit,
    data = mn_training %>% select(-TotalClasses), 
    y = mn_training %>%  pull(TotalClasses),
    label = "rf"
  )

## Preparation of a new explainer is initiated
##   -> model label       :  rf 
##   -> data              :  235  rows  93  cols 
##   -> target variable   :  235  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.2 , task regression (  default  ) 
##   -> predicted values  :  numerical, min =  0.64875 , mean =  4.921705 , max =  43.57892  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -8.43225 , mean =  -0.06638546 , max =  10.42108  
##   A new explainer has been created! 

# compute variable importance
rf_var_imp <- model_parts(rf_explain)

Implications

One note that is critical to keep in mind is that correlation does not imply causation. Although this project looks at connections between various data sets and different variables, we are not suggesting that any of our predictors directly alters computer science course availability. It is possible that that is the case, given the work that we have done, but without an experiment or accounting for potential confounders (other variables that may affect the predictor and outcome variables), we cannot be certain about causation.

Previous work does exist about how disparities in education are related to many of the variables we displayed in our project, such as household income and race. For instance, it has been proven that ACT and standardized test scores show more about family wealth and privilege than actual intelligence or likelihood for success. Therefore, it is logical that the same districts that have high average ACT scores will have high median household incomes due to systemic inequality. Due to the fact that computer science is a newer field, less work has been done specifically about this subject. The rise in available literature on this subject in recent years has also been focused more on college and graduate school, with K-12 education receiving less attention.

That being said, there are many nuances to this issue of course availability and inequality that we could not address within the scope of our project. One variable we looked at was race, and while the connection between race and educational disparities has been studied, that can be difficult to see in some of our work. We hypothesize that there are a few reasons for this. First, Minnesota in general is largely populated by White people. Furthermore, the places with the most racial diversity (near the Twin Cities), are also places with considerable inequality. Without this information, it may seem as though there is correlation between greater diversity, higher median household income, and computer science course availability. However, we cannot make this claim without further investigating how the inequalities within each district play a role.

Along with that, the population size of the districts could affect the outcomes. Districts can encompass many schools, and it is possible that within a district there is variation in demographics and course availability. Future work might include investigating a smaller region to explore some of these nuances in order to better understand the connection between computer science course availability in K-12 education and our other variables.